//============================================================================
// Name        : source_analysis.cpp
// Author      : gabi
// Version     :
// Copyright   : free for all
// Description : Hello World in C, Ansi-style
//============================================================================

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TFile.h"
#include <cmath>
#include "TGraph.h"
#include "TH2F.h"

using namespace std;

int main(int argc, char* argv[]) {
	//TApplication theApp("theApp", 0, argv); //line to allow interactive canvas in the program

	string ifname(argv[1]);
	string ofname = argv[2];
	double maxTheta=atof(argv[3]);
	string ofnameData =argv[4];
	cout << "reading data from " << ifname << endl;
	ifstream ifile(ifname.c_str());
	if (!ifile) {
		cout << "file " << ifname << " not found" << endl;
		exit(1);
	};
	TFile *ofile = new TFile(ofname.c_str(), "RECREATE");
//	if (!ofile) {
//			cout << "file " << ofname << " not opened for writing" << endl;
//			exit(1);
//		};

	string linenames[14] = { "K", "GEN", "Weight", "T(m)", "X(m)", "Y(m)",
			"S(m)", " E(eV)", "Px(eV/c)", "Py(eV/c)", "Ps(eV/c)", "Sx", "Sy",
			"Ss" };
	int nBins[14] = { 10, 10, 10, 10, 1000, 1000, 1000, 1000, 1000, 1000, 1000,
			1000, 10, 1000 };
	//int lowBin[14]={0,0,0,0,0,0,
	double minBin[14];
	double maxBin[14];
	for (int i = 0; i < 14; i++)
		minBin[i] = maxBin[i] = 0;
	double linedata[14];
	int counter;
	counter = 0;
	//for(int i=0;i<14;i++)cout<<linenames[i]<<" ";
	//cout<<endl;
	while (counter < 1000 && ifile) {
		//while (ifile) {
		counter++;
		for (int i = 0; i < 14; i++) {
			ifile >> linedata[i];
			if (linedata[i] > maxBin[i])
				maxBin[i] = linedata[i];
			if (linedata[i] < minBin[i])
				minBin[i] = linedata[i];
		};
//		for (int i = 0; i < 14; i++)
//			cout << linedata[i] << "  ";
//		cout << endl;

	};
	cout << " data lines found:" << counter << endl;
	minBin[7] = minBin[7] / 1000.;
	maxBin[7] = maxBin[7] / 1000.; //E in keV
	for (int i = 0; i < 14; i++) {

		cout << i << " " << linenames[i] << "  " << minBin[i] << "  "
				<< maxBin[i] << endl;
		minBin[i] = minBin[i] * 1.2;
		maxBin[i] = maxBin[i] * 1.2;
	};

	TH1F* hists[14];
	for (int i = 0; i < 14; i++) {
		hists[i] = new TH1F(linenames[i].c_str(), linenames[i].c_str(),
				nBins[i], minBin[i], maxBin[i]);
	};
	TH1F* thetaHist = new TH1F("Theta", "Theta", 10000, 0, 0.2);
	TH1F* thetaSignHist = new TH1F("ThetaSign", "ThetaSign", 100, -0.025, 0.025);
	TH1F* thetaSignHist2 = new TH1F("ThetaSign2", "ThetaSign2", 100, -0.002, 0.002);
		TH1F* thetaSinHist = new TH1F("SinTheta", "SinTheta", 10000, -2, 2);
	TH2F* thetaE = new TH2F("theta_E","theta_E",1000,0,0.25, nBins[7],minBin[7],maxBin[7]);
	TH2F* thetaE2 = new TH2F("theta_E2","theta_E2",1000,0,0.002, nBins[7],minBin[7],maxBin[7]);
	TH2F* thetaE3 = new TH2F("theta_E3","theta_E3",1000,0,0.0004, nBins[7],minBin[7],maxBin[7]);
			double theta;
	theta = 4E-4;
	double thetaArray[40];
	double nThetaArray[40];
	int nThetaArrayC[40];
	double ESum[40];
	double E2Sum[40];
	double Eaverage[40];
	double bw[40];
	for (int i = 0; i < 40; i++) {
		thetaArray[i] = theta / 40 * i;
		nThetaArray[i] = nThetaArrayC[i] = 0;
		ESum[i] = E2Sum[i] = 0;
		bw[0] = Eaverage[i] = 0;
	};

	ofstream ofileData(ofnameData.c_str(),ios::trunc);
	cout<<"writing truncated file to "<<ofnameData<<endl;
if(!ofileData){cout<<"output file for Data not opened"<<ofileData<<endl; exit(1);}
	ifile.clear();
	ifile.seekg(0);

	counter = 0;
	double weight;
	for (int i = 0; i < 14; i++) {
		ifile >> linedata[i];
	};
	weight = linedata[2];

	ifile.clear();
	ifile.seekg(0);
	while (ifile) {
		counter++;
		for (int i = 0; i < 14; i++) {
			ifile >> linedata[i];
			if (i == 7)
				hists[i]->Fill(linedata[i] / 1000., weight); //energy in keV for graphs
			else
				hists[i]->Fill(linedata[i], weight);
		};
		theta = sqrt(linedata[8] * linedata[8] + linedata[9] * linedata[9]);
		theta = atan2(theta, linedata[10]);
		double thetaSign, sinTheta;
		thetaSign = theta * linedata[8] / abs(linedata[8]);
		sinTheta = sin(theta);
		thetaHist->Fill(theta, weight);
		thetaSinHist->Fill(sinTheta, weight);
		thetaSignHist->Fill(thetaSign, weight);
		thetaSignHist2->Fill(thetaSign, weight);
		if(theta<maxTheta){
			//ofileData<<theta<<"  ";
			for(int i=0;i<14;i++)
			ofileData<<linedata[i]<<"  ";
			ofileData<<endl;
		};
		int bin;
		bin = 39;
		//cout<<theta<<"  "<<thetaArray[bin]<<endl;
		while (theta < thetaArray[bin]) {
			ESum[bin] = ESum[bin] + linedata[7];
			E2Sum[bin] = E2Sum[bin] + linedata[7] * linedata[7];
			//cout<<weight<<endl;
			nThetaArray[bin] = nThetaArray[bin] + weight;
			nThetaArrayC[bin]++;
			bin--;
		};
		//cout<<theta<<"  "<<thetaArray[39]<<endl;

		//		for (int i = 0; i < 14; i++)
		//			cout << linedata[i] << "  ";
		//		cout << endl;
		thetaE->Fill(theta,linedata[7]/1000.,weight);
		thetaE2->Fill(theta,linedata[7]/1000.,weight);
		thetaE3->Fill(theta,linedata[7]/1000., weight);
					};
	cout << " data lines histogrammed:" << counter << endl;

	for (int i = 0; i < 14; i++) {
		//hists[i]->Draw();
		hists[i]->Write();
	};
	thetaHist->Write();
	thetaSinHist->Write();
	thetaSignHist->Write();
	thetaSignHist2->Write();
		thetaE ->Write();
	thetaE2->Write();
	thetaE3->Write();

	for (int i = 1; i < 40; i++) {
		Eaverage[i] = ESum[i] / nThetaArrayC[i]/1000.; //keV
		bw[i] = sqrt(
				(E2Sum[i] - ESum[i] * ESum[i] / nThetaArrayC[i])
						/ nThetaArrayC[i]) / 1000.; // /1000 for keV
		bw[i] = bw[i] / Eaverage[i];
	};
	cout << 0 << " " << thetaArray[0] << "  " << nThetaArray[0] << "  " << bw[0]
			<< "  " << Eaverage[0] << endl;
	for (int i = 1; i < 40; i++) {
		//nThetaArray[i] = nThetaArray[i] + nThetaArray[i - 1];
		cout << i << " " << thetaArray[i] << "  " << nThetaArray[i] << "  "
				<< bw[i] << "  " << Eaverage[i] << " " << endl;
	};
	TGraph* gr = new TGraph(40, thetaArray, nThetaArray);
	gr->SetName("flux_angle");
	gr->Write();

	TGraph* gr2 = new TGraph(40, thetaArray, bw);
	gr2->SetName("bw_angle");
	gr2->Write();

	ofile->Close();
}
;
